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1. Motivation and objectives 
The objective of the present research is the development and application of effi- 
cient adaptive numerical algorithms for the study, via direct numerical simulations, 
of active vortex generators. We are using innovative computational schemes to in- 
vestigate flows past complex configurations undergoing arbitrary motions. Some of 
the questions we try to answer are: Can and how may we control the dynamics of 
the wake? What is the importance of body shape and motion in the active control 
of the flow? What is the effect of three-dimensionality in laboratory experiments? 

We are interested not only in coupling our results to ongoing, related experi- 
mental works, but furthermore to develop an extensive database relating the above 
mechanisms to the vortical wake structures with the long-range objective of devel- 
oping feedback control mechanisms. This technology is very important to aircraft, 
ship, automotive, and other industries that require predictive capability for fluid 
mechanical problems. The results would have an impact in high amgle of attack 
aerodynamics and help design ways to improve the efficiency of ships and sub- 
marines (maneuverability, vortex induced vibration, and noise). 

2. Accomplishments 

This is a preliminary report describing our numerical method and presenting 
results of direct numerical simulations of pertinent two-dimensional bluff body flows. 
Our final objective is the simulation of flows past the three dimensional configuration 
shown in Fig. la. At the first stage of this work, we are interested in the respective 
two-dimensioned configuration (Fig. lb). The results of these simulations would 
help us investigate the effects of three dimensionality in experiments and assess the 
role of two-dimensional vortical mechanisms. 

2.1 Mathematical formulation - numerical method 
Two-dimensional incompressible unsteady flow of a viscous fluid may be deter- 
mined by the vorticity transport equation as 

+ u • Yu; = vV 2 u> (1) 

at 

where u(x,f) is the velocity, « = we x = V x u the vorticity, and v denotes 
the kinematic viscosity. For flow around a non-rotating flat plate, translating with 
velocity U b (t), the velocity of the fluid (u) on the surface of the body (x a ) is equal 
to the velocity of the body: 

u(x J} f) = U b (t) 
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Figure 1. Definition sketch. Top view (a), front view (b). 
At infinity we have: 


u(x) -► Uoo as |x| — » oo (lb) 

where Uoo is the free stream velocity. Using the definition of the vorticity and the 
continuity (V • u = 0), it can be shown that u is related to u? by the following 
Poisson equation: 

V 2 u = -V x lj . (2) 

The velocity- vorticity formulation helps in eliminating the pressure from the un- 
knowns of the equations. However, for bounded domains, it introduces additional 
constraints in the kinematics of the flow field and requires the transformation of the 
velocity boundary conditions to vorticity form. 

2.1.1 P article /vortex methods 

The present numerical method is based on the discretization of the above equa- 
tions in a Lagrangian frame using particle (vortex) methods. The vorticity equation 
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Eq. 1 may be expressed in a Lagrangian formulation by solving for the vorticity- 
carrying fluid elements (x a ) based on the following set of equations: 


dx a , j \ 
_ = u(x.,<) 


dw 

dF 


(3) 


^ = t/V 2 


’U) 


In the context of particle methods, it is desirable to replace the right-hand side 
of equations Eq. 3 by integral operators. These operators are discretized using the 
locations of the particles as quadrature points so that ultimately Eq.3 is replaced 
bv a set of O.D.E.’s whose solution is equivalent to the solution of the original s 
of equations. To this effect the velocity field may be determined by the vorticity 
field using the Green’s function formulation for the solution of Poisson s equation 

(Eq. 2). 

u = = / K(x - y) x u>dy + U 0 (x,t) (4) 

27T J 

where U 0 (x,f) contains the contribution from the solid body rotation and U«>, 
and K(z) = z/|z| 2 . The use of the Biot-Savart law to compute the velocity field 

guarantees the enforcement of the boundary condition at infinity. 

The Laplacian operator may be approximated by an integral operator (Degon 
and Mas-Gallic, 1989) as well so that: 


V 2 u> 


j Ge(|x - y|) [w(x) - w(y)] dy 


(5) 


where, in this paper, G t is taken to be G t ( z) = (2/*e>) exp(-\z\ /2c ). The 
boundary condition Eq. la is enforced by formulating the physical mechanism i 
describes. The surface of the plate is the source of vorticity that enters the flow. A 
vorticity flux (dco/dn) may be determined on the boundary in a way that ensures 
Eq. la is satisfied. Here a fractional step algorithm is presented that aUows for the 
calculation of this vorticity flux (see also Koumoutsakos et al. 1993). It shown 
then that this mechanism of vorticity generation can be expressed by an integral 

operator as well: n 

I) - jB^y)^y)dy (6) 

where the kernel H is developed in Section 3. Using the above integral representa- 
tions for the right-hand side of Eq. 3, we obtain the following set of equations. 


dx a 

dt 

du) 

dt 

du 

dt 


/ K(x a - y) x u)dy + U„(x a ,t) 
2n J 

V j Gf(|x a - y|) p(x a ) - w(y)l dy 

J tf(x a , y ) ^(y) ^ 


(V 
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In vortex methods, the vorticity field is considered as a discrete sum of the individual 
vorticity fields of the particles, having core radius e, strength T(<), and an individual 
distribution of vorticity determined by the function so that: 

N 

u>(x,<) = 5^r n (f)r/ e (x -x„(t)) (8) 

When this expression for the vorticity is substituted in Eq. 7, the singular inte- 
gral operators Iv,G are convolved with the smooth function rj ( and are replaced 
by smooth operators K ( ,G f . These integrals are subsequently discretized using a 
quadrature having as quadrature points the locations of the particles. Assuming ! 
that each particle occupies a region of area h 2 and that the shape of the body is ) 
discretized by M panels, then algorithmically the method may be expressed as: 

dxi 1 

~df = A '<( x < ~ x j) + U 0 (x„t) 

j=l 

dr N 

~dt = u ~ r i] G «(l*i-*>l) 

>= 1 

T Tf , duj . 

dt v ^ ^( x, ’ xm ) 3^( x "«) 

m = l 

r,(0) = ui(x,, 0) h 2 I = 1,2,..., N 

The Lagrangian representation of the convective terms avoids many difficulties asso- 
anted with its discretization on an Eulerian mesh such as excess numerical diffusion. 
However the straightforward method of computing the right-hand side of (Eq. 8), 
using (Eq. 9) for every particle, requires 0(N 2 ) operations for N vortex elements! 
Recently fast methods have been developed that have operation counts of 0(N ) 
(Greengard and Rohklin, 1987). In the present scheme the efficient vectorization 
of the O(N) scheme allows for computations with one CPU minute per velocity 
evaluation for one million vortices on a single processor of a CRAY YMP. The ac- 
curacy of the method relies on the accuracy of the quadrature rule as information 
needs to be gathered from the possibly scattered particle positions. The conver- 
gence properties of vortex methods with a finite core dictate that the particles must 
overlap (i.e. h < e) at all times (Beale, 1986). However the local strain field of the 
flow may distort the particle locations, thus producing particle clustering in one 
direction accompanied by an expansion in another direction, similar to that which 
would occur around a hyperbolic stagnation point in a steady flow. When such a 
situation occurs the particle locations have to be re-initialized (remeshed) onto a 
uniform grid while interpolating the old vorticity on the new particle locations. In 
our algorithm we use a nine-point scheme to perform this remeshing, conserving 
the circulation as well as the linear and angular momentum of the vorticity field. 

See Koumoutsakos (1993) for further discussion. 
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2.2 A fractional step algorithm - boundary conditions 
The no-slip boundary condition accounts for the generation of vorticity on the 
surface of the body. The surface of the body acts as a source of vorticity (Lighthill, 
1963), and the task is to relate this vorticity flux on the surface of the body to 
the no-slip condition. This is achieved by appropriately coupling the kinematic and 
dynamic constraints of the flow. The present methodology was first proposed by 
Schmall and Kinney (1974) and is implemented here in a fractional step algorithm. 

In the present formulation Eq. 9 are not integrated simultaneously in time, but 
instead a fractional step algorithm is employed. The governing equations are solved 
via a splitting scheme that accommodates the enforcement of the boundary con- 
ditions. For a non-zero thickness body the enforcement of the no-tangency flow 
implies the enforcement of the no-through flow as well (Koumoutsakos, 1993). This 
is not the case, however, for the zero thickness flat plate, as one has to account 
explicitly for the no-through flow boundary condition. 

Let us assume that at the n-th time step (corresponding to time t — St) the 
vorticity field has been computed (respecting the no-slip boundary condition) and 
we seek to advance the solution to the next time step (time t). The following 
two-step procedure is implemented: 

• Step Is Using as initial conditions /(x) = o?"(x n ,n6f) we solve: 

uj t + u - Vo? = uS7 2 uj in Px[<-ft,<] 

u>(x,t — 6t) = /(x) in V 


where T> denotes the flow domain exterior to the body surface dT>. Particles are 
advanced via the Biot-Savart law, and their strength is modified based on the 
scheme of Particle Strength Exchange. 

Algorithmically then Step 1 may be expressed as: 


dx 
d t 
du^ 
d t 


iz n (x n , n6t ) 
i/V 2 ^ 


(ii) 


At the end of Step 1 a vorticity field (*>[ has been established in the fluid. 

• Step 2: 

The boundary conditions are enforced in this stage by a vorticity (not particle) 
creation algorithm. 

2.2.1 Simply connected domain - flat plate 

Without loss of generality we assume that the plate surface lies along the x-axis 
at y = 0 and between —Lf 2 < x < Lf 2. First, in order to enforce the no- through 
flow boundary condition, a vortex sheet of strength /c(£) is distributed on the plate 
surface. The strength of this vortex sheet is then determined by the solution of the 
following integral equation: 

r jsu = «„(*) 

J — i/2 X ~ £ 


( 12 ) 
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where u„(x) is the velocity normal to the surface of the plate. The above integral 
equation is singular as it admits more than one solution. A unique solution is 
obtained by enforcing the conservation of circulation (note that in airfoil theory k 
is uniquely determined by enforcing the Kutta condition) 



(13) 


Introducing the transformation 6 = cos ~ l ( and using the orthogonality of a 
Fourier series expansion for k( 0), (Batchelor, 1967) over the interval [0, 7r], we obtain 
that at each instant of time: 


k(9) = -(A 0 + 2 sin/?)^£ + + 2 ^ A m sin(mfl) (14) 

m= 1 

where: 

2 f w 

Am = — / u n (8)cos(m6)d0 form = 0, 1, . . . , oo (15) 

n Jo 

and C is determined by enforcing the conservation of circulation (Eq. 13): 

C = § + A, (16) 

A finite number of terms (P) is retained in the infinite series Eq. 14 in the numerical 
implementation of this method. Note also that the expression in Eq. 14 has built 
in the singularity of the vortex sheet at the tips of the plate. The enforcement of 
the no- through flow does not imply the enforcement of the no- slip condition for the 
zero thickness flat plate. A tangential velocity u*(x) may be computed along the 
plate surface due to the vortices in the wake. By taking the limit of Eq. 12 a vortex 
sheet of strength 7 is observed on each side of the plate as: 

j(y = 0~,x) = -^k(x) + ut(x) 

\ (17) 

l(y = 0 + ,a;) = +-/c(z) + u t (x) 

In order to enforce the no-slip condition then, in the context of the present algo- 
rithm we need to eliminate the spurious vortex sheet (k) and the tangential velocity 
on the surface of the plate. 

2.2.2 Doubly connected domain - vortex generators 

For a doubly connected domain as that shown in Fig. lb, the potential flow 
problem needs to be modified. For an open domain such as that exterior to the 
cavity (in the absence of the flat plate), in order to enforce the no- through flow 
boundary condition we require that the tangent flow on the interior surface of the 
boundary vanishes. Unlike the case of a closed body, no additional constraint needs 
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to be imposed for the solution of the respective integral equation (Koumoutsakos 

and Leonard, 1989). . . 

For the combined configuration the resulting integral equation is then: 


,L/2 

J-L/2 


«(0 


L/2 X P 


di + 


r 7ui(o 

J — OO 


ds(x p ) 


-L/2 


- Log\C -x P \dC = u„(x p ) 

(18a) 

f L< 2 duJ ^ p,T ^dx p dT 

l-L/2 On{Xp) 

(18b) 


o l ° 3K -m + = u ' (0 (18c) 

where k(£),x p and 7„,(C),< denote vortex sheets and location of points on the 

surface of the plate and the cavity respectively. 

The above set of equations, when discretized using a panel method, results in 
a well posed system of equations which can be solved iteratively or by direct LU 
decomposition (if storage is not a limiting factor). 

The spurious vortex sheet ( 7 ) that is observed on the surface of the body may 
then be translated to a vorticity flux (Koumoutsakos and Leonard (1994)). 

^(x) = 2(.) on &D (19) 

The computed vorticity flux generates vorticity in the fluid. The vorticity field is 
augmented by this viscous mechanism as described by the following set of equations: 

- vV 2 u2' 2 =0 in x>x[t-6t,t] 

uj' 2 (x,t-6t) = 0 in V ( 20 ) 

du>' 2 _ 7(^1 ) on 

V dn 6t 

Note that the diffusion equation is solved here with homogeneous initial conditions 
as the initial vorticity field was taken into account in the previous substep. The 
solution at Step 2 is a vorticity field lj' 2 , which we superimpose onto the solution of 
Step 1 to obtain the vorticity distribution at the next time step 

w n+1 = +u4 (21) 


3. Results 

We have conducted a computational study of the unsteady flow behind a zero 
thickness plate started impulsively or uniformly accelerated normal to the flow. 
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FIGURE 2. Evolution of the recirculating bubble for Re = 20. : Present 

computations, • ,x : Experiments (Coutanceau and Launay)for blockage ratios of 

0.1 and 0.15 respectively. Theory (Steady State Results): : Dennis and Qiang, 

: Ingham and Tang. 


For the impulsively started plate the present results complement related experi- 
mental works while providing quantities such as the vorticity of the flow field and 
the forces experienced by the body. The development of the flow is similar for all 
the Reynolds numbers that were simulated (20 to 1000). The separating shear layer 
rolls up into a vortex in the lee of the plate, inducing initially a region of secondary 
vorticity. Diffusion acts to increase the width of the shear layer and reduce the 
strength of the vortex resulting in a stable configuration. In Fig. 2 we present the 
results of the present computations for the length of the recirculating bubble for 
Re=20 and compare with existing experimental and theoretical works. The dis- 
crepancies may be attributed to the different treatment of the boundary conditions 
from the calculations of the steady state results and to the finite blockage ratio in 
the experiments. A different behavior is observed for the separating shear layer 
of a uniformly accelerated plate( U = a t). The continuous increase of the shear 
flow overcomes the effects of diffusion, increasing the strength of the separating 
shear layer and inducing a Kelvin- Helmholtz type instability. The wavelength and 
the onset of this instability depends on the acceleration of the plate. The present 
simulations are the first to confirm related experimental evidence on the forma- 
tion of vortex centers along the separating shear layers of an accelerating flat plate 
(Fig. 3). Such undulations have been attributed to experimental defects, but the 
present simulations suggest that this is an intrinsic behavior of the flow. Finally 
the drag coefficient of the plate is shown to scale due to the similarity in the inviscid 
development of the flow. 

An extensive study of these flows for various Re numbers and acceleration rates 
may be found in Koumoutsakos (1994). 
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FIGURE 3. Flow past an accelerating flat plate. Top: Experimental results from 
Lian and Huang (1989) for a - aL 3 /u 2 = 24.51 x 10 6 at T* = at 2 / L = 4.37. 
Bottom: Computational results. Equivorticity Contours for a = 1.68 x 10 6 at 
T* = 6.25. 

4. Future plans 

We are currently at the final stages of development of our code for the study 
of flows involving the doubly connected configuration shown in Fig. lb. We are 
interested in investigating the role of three dimensionality in the experiments and 
determining the role of two-dimensional vortical structures. We will be investigating 
possible control mechanisms with the long-range objective of developing feedback 
control mechanisms. 

Our final goal is full three-dimensional direct numerical simulations to comple- 
ment the ongoing experimental investigations. We may use existing algorithms or 
algorithms under development or proceed to extend the present computational tech- 
nique to three-dimensional flows. In particular, research would be directed in the 
development of fast solutions of integral equations in three dimensions (with pos- 
sible multi-disciplinary applications) and the process of restoring the Lagrangian 
computational grid. Of particular interest would be the coupling of the present 
Lagrangian method with large-eddy simulation (LES) methodologies. 
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